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SUMMARY 

A study was conducted to develop a method for the direct 
integration of a rotor dynanics systen experiencing a blade loss 
induced rotor rub. The approach was first to insure the numerical 
stability of the integration technique: and second, to provide a 
variable time step. This was important because during the rub hi y h 
frequency, vibration components could he excited and smaller time 
steps would he necessary to calculate the vibration components 
accurately. TXirin", other times larger time steps could be used to 
conserve computational time. The method was numerically stable for 
any time step in to a third order integrator. The time step was 
controlled so that the maximum error was less then .01% and the 
probable error was between .001% to .0001 %. 

An ex is tine rotor, (which dynamically simulates a typical snail 
eas turbine) was modeled. The rotor bearing system consisted of a 
shaft with three disks mounted on two preloaied ball hearings, (two 
disi'.s outboard of the bearings) . The bearings were mounted in 
squeeze-filn dampers, which had centering springs. The first three 
critical speeds were calculated to be 7600 , 9203 . and 11200 rpn- all 
three modes are bent-3haft type. °rior to the blade loss simulation 
the rotor was assumed to he balanced and operating at 9 500 rpm . The 
blade loss was simulate:! by an instantaneous application of 5 mils o t 
mass excentricity in the far disk. The rotor rub was simulated by 
surround inq each disk with a shroud that had a 2 mil radial clearance 
and a stiffness of 100,000 ."/in. 

In general both the blade loss and rotor rub phenomenon generate 
hi’h frequency v lbration components. The rotor motion is Initially 
localized but with tine progresses to other Darts of the rotor bv 


means of traveling waves. The traveling waves from several rubs can 
interact with one another causing very complicated rotor motion. Fven 
if there is no rub, (just a blade loss), the traveling wave can cause 
the rotor to beat at a frequency which is the difference between the 
operating and the critical speeds. Rotor rubs generate a frictional 
force vtaich tends to drive the rotor to whirl in a direction opposite 
to the direction of rotation, (backward whirl). For the rotor typical 
of small gas turbines, a snail change in the coefficient of friction, 
(from .1 to .2) caused the rotor to change from forward to backward 
whirl and to theoretically destroy itself in a few rotations. This 
method provides an analytical capability to study the susceptibility 
of rotors to rub induced backward vrtiirl problems. A 10 minute, 
16-raill imeter , color, sound motion picture supplement is available, on 
loan, from the h T ^SA Le'*/is Research Center, that shows the computer 
made motion pictures for the blade loss induced rotor rubs. 


I'rROhTjcTion 

In a typical aircraft gas turbine there are many instances in 
which rotor rubs occur. TWo of the most common are blade tip and seal 
rubs, which are caused by thermal mismatch, rotor imbalance, high "g ,f 
maneuver loads, aerodynamic forces, etc. Current interest in fuel 
efficiency is a consideration which drives the engine design toward 
closer operating clearances. Thus increasing the probablitv of rotor 
rubs. The interaction of a rotor with its case, (rotor rubs) , has 
been studied in ref 1 ami 2. Ref 1 studied a steady state interaction 
between a rotor with a rigid case neglecting friction at the interface 
and Ref 2 studied a steady state interaction between a linear 
flexible rotor and case including friction at the interface. Ref 1 and 
2 did not consider the critical transient situation in which the rotor 
bounces off the case. 

It is known that rotor rubs can have an important effect on the 
rotor dynamics . 'when a rotor rubs on the case , a frictional force is 
generated which can drive a rotor to whirl in a direction opposite to 
the direction of rotation, (backward vhirl) . This frictional force is 
relatively constant up to the backward whirl speed at Uiich the rotor 
mils around the case. Since this rolling contact speed is 
proportional to the rotational speed of the rotor tines the ratio of 
the diameter to the rotor clearance, the v/iirl speei can be hundreds 
of times the rotational speed of the rotor: and thus be potentially 
very dangerous. 

There are two basic methods for studying transient rotor 
dynamics. One is the modal met hoi (ref T and 4) which expands the 
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solution in terras of a few of the lower frequency mode shapes- If the 
transient under study is localized (like a blade loss or a rotor rub), 
the high frequency components are. at least initially, dominant* Thus 
the raodal method is not applicable to this type of transient- The 
other method involves the direct integration of the equations of 
motion, which can be done in either of two ways, explicit or implicit 
integration. For example, ref 5 used explicit integration of the 
equation of motion, but this solution is plagued with numerical 
stability problems- Further, ref 6 showed that explicit integration of 
the equation of motion was unstable when the product of the critical 
frequency (for any mode numerically possible) and the time step was 
large- Therefore, the explicit integration can only be done for simple 
rotors - 

In contrast, the implicit integration tends to be stable (ref 7 
and 8) but it requires the solution of a large number of nonlinear 
simultaneous equations at each tine step. 'He f 9 used a technique 
similar to ref 7 except that it was applied directly to the second 
order equation of motion- Ref 9 also noted that the generalized 
forces on a rotor were functions of the generalized position and 
velocity of the point where the forces were applied and its nearest 
axial neighbors* Tills allowed the variables to be arranged so that the 
Jacobian of the set of nonlinear equations was block tridiagonal * 
Therefore, computing tine became proportional to the number of 
elements in the rotor dynamics node! rather than to the cube of the 
number of elements. The objective of this study is to refine the 
method used in ref 9 to include an automatic time step routine and 
then apply the technique to study blade loss induced rotu^ rubs* The 
automatic time step routine is necessary so that the time step can be 
varied as the rotor impacts the case- Also, the numerical stability 
of the method used in ref 9 will be investigated. 


SYMBOLS 


a reference amplitude 

c radial clearance 

E absolute error estimate 

F force 

0 order of error in Taylor series 

q order of Taylor scries 

r radial displacement 

S stability matrix 

t t ime 

At time step 

u def ined in eq • ( 4 ) 


z Independent variable 

a given set of constants 

C damping ratio 

X eigenvalue of stability matrix 

p coefficient of friction 

m frequency 


ANALYSIS 


N^mericj^^ i nteg ration * 

Given an arbitrary vector function ^(t) whose derivatives exist 
vt)» a Taylor series expansion can be written: 


q-k 


( 1 ) 


j-0 

■+ 

with remainder of order If the arbitrary function is chosen as 


? - iMlltOO 

k ak! 

the Taylor series for this function becomes 

q 


J k (t ♦ At) - £ + 3 q 


j-0 

where the binomial coefficients are defined as 




for ] ik 


/j\.Jk:(j - k)! 

' / l 0 for j < k 


( 2 ) 


(3A) 


( 3H) 


If the form of the remainder is chosen as 


0 - a, u 

q k 


the Taylor series becomes 


Zj,(t + Lt) 
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J-0 


(t) + o^u 


(4) 
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vhere the alphas are given in ref 7 and u can be determined from the 
equations of motion at the advanced time- The form of the aet of the 
equations of motion at the advanced time is: 

£F(r, r, r, t ■+ At) - 0 < 6 > 
From the definition of z . the various derivatives become: 


Substitutin'’ for the various derivatives into the equations of motion; 
“ values at the previous tine, result in the epuations 

of motion being a function of: 

£f(u, t + At) - 0 (8) 

This set of equations can be solved for u and. fron this value of 
the remainder can be used as an error estimate to control the time 
step. From the definition z x represents a nond imensional form of ry 
Therefore an estimate of the maximum absolute error is : 


where I fi,| | . the vector norm is the maximum component of u- Toe 
computer code used in ref 9 was modified to include the followin 3 
automatic time ster althroqim- If E>-01V re-do the calculation with 
the time steo reduced by a factor of 10- If .01Z>E>.001S. accept the 
calculation but decrease the time step by a factor Oi L ■ i 
001 1>T>. 0001 S. accept the calculation and maintain t e sane me 
step If .0001%>r. accept the calculation but increase the time step 

bv a factor of 2- 
Nuner ic al stabilit jy : 

The analysis of the stability of the numerical integration 
technique assumes a model of a rotor bearing system that is l*"®** 1 "* 
at some instant of time- The homogeneous equation of motion for an. 

mode is 

r + 2u)^r + w 2 r - 0 U°) 

where omega is the natural frequenev and zeta is the damping lo for 
the mode. For every mode that is numerically possible, with 


nonnegative damping ratio, the amplitude oust cither remain constant 
or decay in time- The numerical integration is defined as unstable if 
the amplitude grows in time. 

Substituting the Taylor series for the various derivatives into 
the modal equation of motion at the advanced time results in; 
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V' pd - 1) + 2 iu.' At fe 
/ ,« [2a 2 + 20^ At ;+ a Q (w A 
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For this value of u, the Taylor series expresses the solution at the 
advanced time in terms of the solution at the present time as: 


Z k (t + At) 


XffiP 


j-0 


j (j - 1) + 2ju At t + (uj At) 


2a 2 + 2a x u At X, + QqCu At)' 



Z . (t) (12) 


Defining the matrix S to be: 


/j\ a k[ j(j ~ X) 4 2ju Lt 1+ (ai A ° 2 ] 

S kj \ k ] j2a 2 + 2^ At ;+ a Q (w At) 2 


(13) 


and the vector Z whose kth element is z k - results in the finite 
difference equation : 


l(t + At) - Sl(t) 

This equation has a solution of the form; 

l(t + At) - Xl(t) 
where lambda is an eigenvalue of : 

S t - \t 


(14) 


(15) 


(16) 


If the |X | >1, the amplitude grows and the method is numerically 
unstable • 


Rub m odel : 

The interaction of a rotor with its case is a complicated 
phenomenon. It can involve non-linear deformation of both the rotor 
and the case. Rotor-case rubs were experimentally studied in ref 10. 
Analytically only simple rotor-case rub models are available; 
therefore, the case was assumed to be linear with dry friction 
interaction with the rotor* The radial and tangential forces on the 
rotor are then: 


F r • °. F 9 - 0 

! r | < C 

( 1 7A ) 

F r - -k(|I| - C), F 0 - 

|r| > C 

(17b) 


RESULTS ALT) DISCUSSION 

The numerical method of ref 9 employed a second order integrator 
with a constant time step. However, to study blade loss induced rotor 
rubs, it is necessary to modify the nethod of ref 9 to include higher 
order integrators with an automatic time step routine- The autonatic 
time step routine is necessary so that the tine step can be varied as 
the rotor impacts the case- In order to calculate high frequency 
components accurately, the tine step must be less than the period of 
the high frequency component- Vhen only low frequency components are 
important the time step can be increased to decrease computing time- 
Th ^ algorithm discribed in the analysis section keeps the maximum 
error in the displacement at less than .OKI- It tries to maintain the 
error between -001% to -COOl', by either decreasing or increasing the 
time step. 

Another way to decrease computing time is to use a higher order 
integrator. Ref 7 studied the numerical stability of up to a sixth 
order integrator applied to a first order differential equation- The 
numerical stability of these integrators applied to a second order 
differential equation was given in the analysis section. The 
numerical stability of an integrator is based on modal rotor dynamics 
analysis. If the integrator is applied to a node which is not driven 
and has damping, the amplitude must decay in tine. Figure 1 shows a 
stability map for the integrators used in ref 7 applied to a second 
or der differential equation. The abscissa is the damping ratio and the 
ordinate is the product of the time step and natural frequency for the 
mode. The stability map has contours on it for which the amplitude 
does not change from one time to the next- On one side of the contour 


the amplitude grows (unstable region), and on the other side it 
decays, (stable region). 
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Figure 1 shows the stable regions for a fourth through 6ixth 
order Integrator. The second and third order Integrators were stable 
everywhere. For the regions irtiere the Integrators were unstable, the 
anplitude grew by a few percent per time step. It would take on the 
order of a hundred time steps for the amplitude to double, and it 
would take on the order of a thousand time steps for the anplitude to 
increase by a factor of a thousand. Due to round off errors, every 
mode that is numerically possible in the rotor dynamics model, has a 
finite amplitude. These amplitudes may be small- but if they are in an 
unstable region, in a few thousand time steps they can become very 
large. For this reason, only the second and third order integrators 
were used. This is still a vast improvement over other types of 
integrators such as the one used by NASTRAN. NASTRAN uses an implicit 
form of the Nevmark-Beta integrator, ref 8. This integrator is second 
order and does not have an error estimate. 

The rotor-bearing system described in ref 11. (which dynamically 
simulates atypical small gas turbine), was used as the example 
problem. This rotor bearing system consisted of a shaft with three 

disks mounted on two axially preloaded ball bearings (fig 2). In this 
rotor-bearing system the bearings were mounted in squeeze-film damper 
journals, and the Journals had centering springs. 

The first three critical speeds for the rotor bearing system 
without oil in the dampers are shown In figure 3- Note that all the 
modes are bent- shaft modes. The ’•classical" hierarchy only applies 
to stiff shafts- therefore, the classical node shapes do not 
characterize the actual mode shapes. The first mode, about 7600 rpm. 
classically would be the cylindrical node. But in this case, it has a 
large amount of bending outward near the shaft center. The second 
mode, about 9200 rpm, classically would be the conical mode. In this 
case, it has a slight amount of bending outward near the shaft ends. 

he third mode, about 11200 rpn. classically would be the bending 
mode. In this case, it has a large amount of bending throughout the 


The rotor-bearing system was modeled by using 23 elements. Prior 
to the blade loss simulation the rotor was assumed to be balanced and 
operating at 9500 rpn. The blade loss was simulated by an 
Instantaneous application of 5 mils of mass excentricity in the far 
disk. The equations of motion for this system were directly 
Integrated by the method used in ref 9 with a variable time step. The 
output was interpolated to equal time steps (100 time steps per shaft 
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rotation), and displayed on a CRT. figure 4. The display showed an 
oblique view of the rotor bearing system, with the bearing center line 
as the oblique axis. The transverse vibration is indicated by the 
position of the rotor centerline. The scale of the transverse 
vibration exaggerates the amplitude of the vibration. The display on 
the CRT was photographed at each time step. These photographs were 
then shown as a motion picture. 

Figure 5 shows the superposition of the first ten frames of the 
blade loss simulation without a rub. Initially the rotor, the 

bearing, and the mass center line coincided. After the blade loss, a 
traveling wave starts at the blade loss disk and travels down the 
rotor. During the time high frequency components are dominant, 

(because the rotor as a whole is not moving)- A model analysis which 
only uses the lower modes cannot discribe the notion during this time 
period . 

Figure 6 shows the position of the rotor for the first six 

rotations of the rotor after blade loss without a rub taking place- 
During the first rotation, the blade loss disk spirals out. During 

the second rotation, the disk on the other end of the shaft spirals 

out- During the third rotation, the center disk spirals out- After 
this the envelope of the rotor positions, seems to oscillate in a 
conical fashion, with a frequency of about 1/4 operating speed- This 
beating seens to be at a frequency difference between the operating 
speed and the 1st critical speed. (Ref 12 experimentally showed a 
similiar beat frequency betwen the operating speed and the critical 
speed.) During this time the rotor shape resembles the third critical, 
except that the bearing center line is not in the plane of the rotor- 
The maximum amplitude occurred on the blade loss disk on the sixth 
rotation and on the opposite disk on the fourth rotation. The 
conclusion drawn from this figure is that if there is clearance space 
down the rotor and a rub occurs, it does not necessarily occur at the 
blade loss disk first. 

The rotor-case rub was simulated by surrounding each disk with a 
shroud that had a 2 vail radial clearance and a stiffness of 
100, 000/- /in . The rub was induced by a repeat of the blade loss 
simulation with the clearance restrain. Two rub simulations were run. 
one with, a coefficient of friction of .1 and the other with a 
coefficient of friction of .2. 

Figure 7 shows the first 6 rotations of the shaft after blade 
loss for a coefficient of friction of .1. During the first shaft 
rotation the blade loss disk spirals outward and bounces off the case 
four tines. Each collision of the rotor with the case sends out 


traveling wves down the rotor. These waves Interact with each other 
causing the envelope of the rotor motion to be very complicated. Ch 
the second shaft rotation both outboard disks bomce off the case four 
times. As the rotor continues to turn the orbit becomes more 
circular. That is, the rotor-case interaction becomes less of a 
bouncing nature and more of a continuous contact. The envelope of 
the rotor motion seems to be oscillating in a conical nature; but both 
outboard disks see.o to remain in contact with the case. Tie rotor 
continues to Oiirl about the bearing centerline in the rot. .ional 
direction (forward whirl) . The frictional drag forces are not lane 
enough to drive the rotor into backward whirl. 

Figure 3 shows the first 4 rotations of the shaft after blade 
loss for a coefficient of friction of .2. The motion of the rotor on 
the first rotation is similar to the .1 coefficient of friction case. 
On the second rotation, the blade loss disk has a very hard collision 
with the well, causing the rotor to bend considerably. On the third 
rotation the rotor whirl direction changes from forward to backward 
whirl and the rotor whirl begins to accelerate in the backward 
direction. On the fourth rotation . the rotor motion becomes very 
iarge and it continues to grow on succeeding rotations. 


This example problem has shown that snail 
coefficient of friction, (from . 1 to .2) can change 
to a blade loss condition from a relatively safe 
catastrophic response. R> r seal rdbs the coefficient 

probably between .1 to .2. For blade tip rubs, this 
accurate. This type of rih involves material renewal 
and cr non-elastic defoliations. If this model were 
general manner, then the coefficient of friction wo 
greater then .2. 


changes in the 
a rotor response 
response to a 
of f ric tion is 
rub nod el is no t 
, phase charts ft-, * 
to be useJ in a 
uld probably be 


In conclusion, this computer code allows us 
induced rotor rii^s and displays the rotor motion 
fomat. A 10-ninute, lb-mill imete r , color, 
supplement is available, on loan, that shows the 
picture for the blade loss induced rotor rdbs. 


to look at blade* loss 
in a mot ion- pic tu re 
so und notion- pic ture 
computer male notion 


RIMMAS Y OF RESULTS ANh CONCLUSIONS 


A method for direct integration of a rotor dynamics 
experiencing a blade loss indeed rotor rub was develop 
following conclusions were drawn: 

1. The method whs nunerlcally stable for any ti-e step U n to 
order integrator. r 
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2- The time step was controlled so that the raaximun error was less 
then .012 and the probable error was between .001% to .0001%. 

3- For the rotor typical of small pas turbines a snail change in the 
coefficient of friction, (from .1 to .2). caused the rotor to change 


from forward 

to backward 

whirl and 

to destroy Itself 

in a 

few 

rotations . 
This method 

provides an 

analytical 

capability to 

study 

the 


susceptibility of rotors to rub induced backward whirl problems- 
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Figure 4 - Oblique view of rotor centerline whirling 
about the bearing centerline. 



Figure 5. - Initial movement of rotor after 
blade loss. 










